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We study absorbing-state phase transitions in two-dimensional Voronoi-Delaunay (VD) random 
lattices with quenched coordination disorder. Quenched randomness usually changes the criticality 
and destroys discontinuous transitions in low-dimensional nonequilibrium systems. We performed 
extensive simulations of the Ziff-Gulari-Barshad (ZGB) model, and verified that the VD disorder 
does not change the nature of its discontinuous transition. Our results corroborate recent findings 
of Barghatti and Vojta [Phys. Rev. Lett. 113, 120602 (2014)] stating the irrelevance of topological 
disorder in a class of random lattices that includes VD and raise the interesting possibility that dis¬ 
order in nonequilibrium APT may, under certain conditions, be irrelevant for the phase coexistence. 

We also verify that the VD disorder is irrelevant for the critical behavior of models belonging to the 
directed percolation and Manna universality classes. 

PACS numbers: 02.50.Ey, 05.70.Ln, 64.60.Ht, 75.40.Mg 


I. INTRODUCTION 

Nonequilbrium phase transitions from an active (fluc¬ 
tuating) to an inactive (absorbing) phase in spatially ex¬ 
tended systems is a topic of broad interest [TH5] . The so- 
called absorbing-state phase transitions (APTs) arise in a 
wide variety of problems as, for example, heterogeneous 
catalysis |1], interface growth [^, population dynamics 
and epidemiology [5]. Recent experimental realizations 
in turbulent liquid crystals , driven suspensions [8] and 
superconducting vortices [S] highlight the importance of 
this kind of transition. 

In analogy with equilibrium phase transitions, it is ex¬ 
pected that continuous APTs can be classified in univer¬ 
sality classes mm- Generically, single-component sys¬ 
tems with short-range interactions exhibiting a contin¬ 
uous APT, in the absence of extra symmetries or con¬ 
servation laws, belong to the directed percolation (DP) 
universality class [iniiii], but other robust classes emerge 
when multiple absorbing states and conservation laws are 

included mm- 

Of particular interest is how spatially quenched dis¬ 
order affects the critical behavior of an APT. In real 
systems, quenched disorder appears in the form of im¬ 
purities and defects [T^]. On a regular lattice, quenched 
disorder can be added in the forms of random deletion 
of sites or bonds HMH] or of random spatial variation 
of the control parameter umn- In all the cases above, 
quenched randomness produces rare regions which are lo¬ 
cally supercritical even when the whole system is subcrit- 
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ical. The lifetime of active rare regions is exponentially 
long in the domain size. The convolution of rare region 
and exceedingly large lifetimes can lead to a slow dy¬ 
namics, with non universal exponents, for some interval 
of the control parameter < A < Ac where Ac°^ and 
Ac are the critical points of the clean and disordered sys¬ 
tems, respectively. This interval of singularities is called 
Griffiths phase (GP) [22]. This GP behavior was verified 
in DP models with uncorrelated disorder irrespective to 
the disorder strength and corresponds to the universality 
class of the random transverse Ising model [MUMISI- 

These findings are in agreement with the heuristic 
Harris’ criterion [27], which states that uncorrelated 
quenched disorder is a relevant perturbation if 

di^± < 2 , ( 1 ) 

where d is the dimensionality and h'± is the correlation 
length exponent of the clean model. Note that in DP 
this inequality is satisfied for all dimensions d < 4, since 
i/j. = 1.096854(4), 0.734(4) and 0.581(5), for d = 1, 2 and 
3, respectively P5H5D] . In the opposite way, simulations 
of the continuous APT in models with a conserved field 
in the Manna universality class |31j , considering uncorre¬ 
lated lattice dilution below the lattice percolation thresh¬ 
old, provide strong evidences that this kind of disorder 
is irrelevant although the Harris criterion is satisfied for 
d < 4 IM34]. 

For equilibrium discontinuous phase transitions the 
Imry-Ma criterion [33 [311 governs the stability of macro¬ 
scopic phase coexistence and disorder destroys phase co¬ 
existence by domain formation in dimensions d < 2. If 
the distinct phases are related by a continuous symme¬ 
try the marginal dimension is d = 4 |36j . Therefore, 
first-order phase transitions become rounded in presence 
of disorder for d < 2. 

Recent numerical results provide evidences that the 
Imry-Ma argument for equilibrium systems can be ex- 
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tended to non-equilibrium APTs: Irrespective to the un¬ 
correlated disorder strength, Buendia and Rikvold [371- 
I39j reported that the absorbing discontinuous transition 
in the Ziff-Gulari-Barshad (ZGB) model for heteroge¬ 
neous catalysis turns to a continuous one (see also the 
discussion in [40]). Analogous behavior was observed 
more recently by Martin et al. m for a two-dimensional 
quadratic contact-process [H]. 

Another important question is the role played by dis¬ 
order inherent to the underlying connectivity in a non¬ 
periodic, random structure of integer dimension as the 
random lattice generated by the Voronoi-Delaunay (VD) 
triangulation [43] . This random lattice can be gener¬ 
ated from a random (uniform) distribution of N points 
in a unitary square region. The triplets that can be cir¬ 
cumscribed by a circle that does not overlap with any 
other point form a triangulation. The result is a two- 
dimensional connected graph with a Poissonian distri¬ 
bution of connectivity with average degree q = 6 [H]. 
This lattice plays an important role in the description 
of idealized statistical geometries such as planar cellular 
structures, soap throats, etc. [151 Hi] . 

Recently, it was found that such a kind of VD disor¬ 
der does not alter the character of the APT exhibited 
by the clean contact process (CP) [45], which is a proto¬ 
typical model in the DP universality class. These results 
are in evident contrast with those for uncorrelated disor¬ 
der which leads to an infinite-randomness critical point 
and strong GPs [l8l [26] . In order to determine the rel¬ 
evance of the disorder in these cases, we can apply the 
heuristic Harris-Luck criterion |46] . in which the regular 
critical behavior remains unchanged when the wandering 
exponeniQw does not exceed a threshold value given by 

10 ^ = 1 -^. ( 2 ) 

dvi_ 

For independent dilution, uj = 1/2, and Luck’s expression 
reduces to the Harris criterion. 

Former numerical estimates of wandering exponents 
for VD triangulations indicated a value close to indepen¬ 
dent dilution w = 1/2 [47]. So, the clean critical behavior 
observed for CP on VD lattices posed doubts on the va¬ 
lidity of Harris-Luck criterion for DP class [JSj. This 
inconsistency was recently unfolded [35] with the deter¬ 
mination of the correct wandering exponent of VD lat¬ 
tices as w = 1/4 in d = 2 implying a criterion > 2/3, 
not i'± > 1, for a clean critical behavior. 

In the present work, we investigate the role played by 
the disorder of VD lattice on the phase coexistence of 
ZGB model. We provide evidences that the VD topolog¬ 
ical disorder does not destroy the phase coexistence and 
thus permit discontinuous phase transitions. We comple¬ 
ment the paper with more evidences for the irrelevance of 


^ The wandering exponent is associated to the decay of deviations 
from the average as a function of patch sizes where the averages 
are computed. 


VD disorder for continuous APTs belonging to the DP [5] 
and Manna [ST] [49] universality classes. 

The reminder of this paper is organized as follows. In 
the next section, we review the models definitions and 
details of the simulation methods we used. In Sec. HI, 
we present our results and discussions. Sec. IV is devoted 
to summarize our conclusions. 


II. MODELS AND METHODS 

We constructed the Voronoi-Delaunay lattice with pe¬ 
riodic boundary conditions, following the method de¬ 
scribed in [^. For sake of simplicity, the length of the 
domain where N node are randomly distributed will be 
expressed in terms of L = '/N. 


A. Discontinuous APT: 

The ZGB model [3], a lattice gas model introduced to 
investigate the reaction of CO oxidation on a catalytic 
substrate, follows the Langmuir-Hinshelwood mecha¬ 
nism, 

COgas + * —>■ COads 

02g as -|- 2 * —>■ 20ads 

COods + Oads CO 2 -|- 2*, 

where * denotes an empty site, and subscripts indicate 
the state (gaseous or adsorbed) of each species. The 
02 gas dissociates at surface, and requires two empty sites 
to adsorb, while CO requires only one site to adsorb (the 
model is also called the monomer-dimer model). The 
product CO 2 desorbs immediately on formation. COgas 
molecules arrive at rate Y per site while O 2 arrives at rate 
(1 — Y), with 0 < V <1. Varying the control parame¬ 
ter V, the model exhibits phase transitions between an 
active steady state and one of the two absorbing or “poi¬ 
soned” states, in which the surface is saturated either by 
oxygen (O) or by CO. The first transition (0-poisoned) is 
found to be continuous while the second (CO-poisoned) 
is strongly discontinuous. 

The computer implementation is the following: With 
a probability V a CO adsorption attempt takes place, 
and with a complementary probability (1 — V) an O 2 
adsorption attempt takes place. In the former case, one 
site is randomly chosen. If the site is occupied, either by 
O or CO, the attempt fails. If it is empty but one of its 
first-neighbors is occupied by an O, both sites become 
empty (O and CO react instantaneously). Otherwise, 
the site becomes occupied by an adsorbed CO molecule. 
Analogous procedure is followed for an O 2 adsorption 
attempt, but in this case we have to choose at random 
a pair of first-neighbors sites, and check for the opposite 
species in all remaining nearest-neighbors of the target 
pair. 
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B. Continuous APT: 


The CP [5T] is the prototypical model of the DP class, 
and is defined on a lattice with each site either active 
(tJi = 1) or inactive {ai = 0). Transitions from active 
to inactive occur spontaneously at a rate of unity. The 
transition from inactive to active occurs at rate ajX/kj, 
for each edge between active nearest neighbors j of site 
i. The computer implementation of CP in graphs with 
arbitrary connectivity is as follows mm- An occupied 
site is chosen at random. With probability p = 1/(1 + A) 
the chosen particle is removed. With the complementary 
probability 1 — p = A/(l + A), a nearest neighbor of the 
selected particle is randomly chosen and, if empty, is oc¬ 
cupied, otherwise nothing happens and simulations runs 
to the next step. Time is incremented by St = 1/n, where 
n is the number of particles. So, the creation mechanism 
in CP effectively compensates the local connectivity vari¬ 
ation with a reduction of the spreading rate through a 
particular edge inversely proportional to the connectiv¬ 
ity of the site that transmits a new particle. If we modify 
these rules to create offspring in all empty nearest neigh¬ 
bors of the randomly chosen occupied site we obtain the 
A model [52155] (in the A-model occupied sites becomes 
empty at unitary rate, as in the CP). This means that 
sites with higher coordination number produce more ac¬ 
tivity when compared with CP, enhancing possible “rare 
region effect” [22 [54] . Since contagion occurs more read¬ 
ily in the A model than in the CP, the critical creation 
rate Ac is smaller but the two models share the same 
critical behavior of the DP universality class [55] . 

The Manna model |3I|, a prototypical model of the 
Manna class and introduced to investigate the dynamic 
of sandpiles in the context of self-organized criticality, is 
defined on a lattice where each site assumes integer val¬ 
ues (mimicking the number of “sand grains” deposited 
on the substrate). In the version we investigate an un¬ 
limited number of particles per site is permitted. Sites 
with a number bellow a threshold A/, = 2 are inactive 
while those where this number is equal to or larger than 
Nc are active. The active sites redistribute their parti¬ 
cles among its nearest neighbors chosen at random, gen¬ 
erating a dynamics that conserves the number of parti¬ 
cles when considering periodic boundary condition. The 
Manna model exhibits a continuous phase transition from 
an active to an inactive state depending on the control 
parameter p that is given by the density of particles on 
the lattice [551 • The absorbing stationary state, where 
all sites have a number of particles bellow Nc is char¬ 
acterized by an infinite number of configurations. The 
computer implementation is analogous to that of CP: one 
active site i {Ni ^ Nc) is randomly chosen. Each of the 
Ni particles is sent to a randomly chosen nearest neigh¬ 
bor irrespective of its state. The site i becomes empty 
(inactive) and the nearest neighbors with Nj = Nc — 1 
particles that received a new one are activated. 


C. Simulation methods 

The central method we used involves the quasi¬ 
stationary state, in which averages are restricted to sam¬ 
ples that did not visit an absorbing state [Tj. To per¬ 
form the QS analysis we applied the simulation method 
of Ref. [57] . The method is based in maintaining, and 
gradually updating, a set of configurations visited during 
the evolution; when a transition to the absorbing state 
is imminent, the system is instead placed in one of the 
saved configurations. Otherwise the evolution is exactly 
that of a conventional simulation [58] . Each realization of 
the process is initialized in an active state, and runs for at 
least 10® Monte Carlo time steps. Averages are taken in 
the QS regime, after discarding an initial transient of 10^ 
time steps or more. This procedure is repeated for each 
realization of disorder. The number of disorder realiza¬ 
tion ranged from 20 (for the largest size used, L = 2048) 
to 10®. Another important quantity is the lifetime in the 
QS regime, r. In QS simulations we take r as the mean 
time between successive attempts to visit the absorbing 
state. 

For discontinuous APTs, we estimated the transition 
point through the jump in the order parameter and the 
finite-size scaling of the maximum of the susceptibility. 
In DP class the spreading analysis starting from a sin¬ 
gle active site (a pre-absorbing configuration) is very ac¬ 
curate and computationally efficient method [1]. For 
Manna class, spreading analysis is more cumbersome [2 
due to infinitely many pre-absorbing configurations. So, 
we proceeded using dimensionless moment ratios analy¬ 
sis in the QS state [53] . which are size-independent at 
criticality. Here, we analyze the critical moment ratio 
m = (p^)/(p)^, which assumes a universal value me at 
the clean critical point. 

III. RESULTS 
A. Discontinuous APT 

First order transitions are characterized by a disconti¬ 
nuity in the order parameter and thermodynamic densi¬ 
ties, with an associated delta-peak behavior in the sus¬ 
ceptibility [2. However, at finite volume thermodynamic 
quantities become continuous and rounded. According 
to the finite-size theory, rounding and shifting of the co¬ 
existence point scale inversely proportional to the sys¬ 
tem volume L'^ [52. Although there is no established 
similar scaling theory for nonequilibrium systems yet, 
some studies show evidences of an analogous behavior 
for APTs [5nH52 . 

Quasi-stationary analysis remains useful in the context 
of discontinuous APTs [52. Considering the QS simula¬ 
tions we observe a discontinuous phase transition from a 
low-density to a poisoned (absorbing) CO state, as shown 
in Fig. [^instead of a rounded (continuous) transition ex¬ 
pected for APTs in the presence of relevant disorder m- 
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The inset of Fig. shows the QS probability distribution 
for the density of active sites near the transition. We 
clearly observe a bimodal distribution, which is a hall¬ 
mark of discontinuous phase transition |62j . 



Y 


FIG. 1. (Color online) QS Density of CO sites in ZGB 
model on triangular and Voronoi lattices of linear system size 
L = 100, showing a discontinuous phase transition to the CO- 
poisoned absorbing state close to F = 0.56. Inset shows the 
QS distribution for Y — 0.5590 where the APT takes place. 

In analogy to equilibrium first-order phase transition, 
where at the transition point a thermodynamical poten¬ 
tial (such as the free energy) is equal for both phases [59], 
we can define the coexistence value of the order param¬ 
eter in which the area under the peaks of the QS dis¬ 
tribution related to each phase (active and absorbing) 
are equal [62]. The intercept of the linear fit from this 
equal histogram method yields 1). = 0.55928(3). Such a 
value is very close to the coexistence value W = 0.5596(5) 
we found for the regular triangular lattice (see inset of 
Fig.[|). 

The location of the maximum of the susceptibility x^ 
defined as variance of the order parameter y = L‘^((p^) — 
(p)^), scales as in a discontinuous APT Fig¬ 

ure]^ shows the finite-size scaling of the transition point 
which clearly scales inversely to the volume, confirming 
again the disordered lattice does not alter the discontin¬ 
uous nature of the transition. 

Further evidence of discontinuity of the phase tran¬ 
sition in presence of quenched coordination disorder, is 
shown in Fig. [^ Using conventional simulations, we ob¬ 
serve the system bistability around the transition point: 
depending on the initial density, a homogeneous steady 
state may converge either to a stationary active state of 
high CO 2 production (and small CO density) or to the 
CO-poisoned (absorbing) state. These results contrasts 
with those for uncorrelated disorder, for which no matter 
its strength, the discontinuous transition is replaced by 
a continuous one. 



FIG. 2. (Color online) Quasi-stationary susceptibility on the 
ZGB model on VD lattices as a function of Y for different 
sizes. Inset: finite size scaling for the susceptibility maxima 
in the range L = 40 to 320. 



FIG. 3. Density of CO as a function of time for distinct 
initial conditions close to the transition point. Initial densities 
Pco = 0.0, 0.1, 0.2, ■ ,0.9, from bottom to top. Linear 

system size L =100 and Y = 0.5560. 


B. Continuous APT 

The spreading analysis for the A model on VD lat¬ 
tices using the mean number of active sites against time, 
with a single occupied site as initial condition, provides 
a critical value Ac = 0.322430(5), which is smaller than 
Ac = 0.34047(1) found for the regular triangular lattice 
with (7 = 6 . This difference is more significant than that 
obtained for CP for these same lattices [IS] showing that 
effect of disorder in A model is stronger than in CP. 
However, the critical behavior remains that of the clean 
system, exhibiting clear power laws with spreading ex¬ 
ponents very consistent with the DP class (results not 
shown). 

Figure]^ shows that at the critical point we found the 
QS density p decays as a power law, p ~ with 
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/3/i/_L = 0.79(1). Besides, we observe that the lifetime 
of the QS state also follows a power-law at criticality, 
with T ~ L^, z = 1.73(5). Both values of the exponents 
are close to the DP ones of /3 /j^_l = 0.797(3) and z = 
1.7674(6) [2]. The inset of Fig. shows the ratio m = 
{p^)/{p) around the criticality for varying system sizes. 
From these data we found rric = 1.33(1), in agreement 
with the value rric = 1.3257(5) found for DP class in two 
dimensions [53| . All results presented here confirm the 
irrelevance of disorder of the VD lattice for the critical 
behavior of the A model. 



FIG. 4. (Color online) FSS of the critical A model. Main: 
Quasistationary density of active sites p (stars) and lifetime 
of the QS state r (crosses) as a function of the system sizes 
L for A = 0.32243. Inset: Quasistationary moment ratio m 
versus 1/L, for A = 0.32238, A = 0.32242, A = 0.32246, from 
top to bottom. 


Lets now turn our attention to the Manna class. The 
correlation length exponent = 0.799 [S3], larger than 
the DP value 0.7333, makes the modified Harris-Luck 
criterion modified criterion {d+l)i'± < 2 still not fulfilled 
for VD lattices [IBj. Critical point determination using 
moment ratios is shown in the inset of Fig. resulting 
in the estimate Pc = 0.688808(2) that is smaller than 
the triangular lattice threshold Pc = 0.69375(5). The 
critical moment ratio is rric = 1.35(l)jwhich agrees with 
the value we found for square lattice^ruc = 1.348(7) at 
the threshold Pc = 0.716957(2). The critical exponents 
we obtained using L ^ 256 were = 0.78(1) and 

= 1.54(2) are also in striking agreement with the 
Manna class exponents = 0.80(2) and = 

1.53(5). 

It is known that critical exponents and moment ratios 
of the Manna class in d = 2 obtained via QS analysis 


^ Our estimate of m does not agree with that of Ref. m where 
a restricted version of the Manna model, in which > 2 is 
forbidden, was considered. 



FIG. 5. (Color online) Critical Manna model on VD 
lattices. Main: Critical density of active sites (crosses) 
and lifetime (stars) against lattice size. Inset: Moment 
ratio m = (p^)/(p)^ against inverse of size for p = 

0.688800, 0.688805, 0.688810, and 0.688815 from top to bot¬ 
tom. 


are hardly distinguishable form DP class HISS]. In or¬ 
der to provide a more incisive verihcation that Manna 
model on VD lattice has exponents different from DP we 
considered density around the critical point, which scales 
as [55] 


where A = p — pc- This implies that 


d\np 


dv 




(3) 

(4) 


can be used to obtain the exponent explicitly. 
Similarly, for the moment ratio we have m{A,L) = 
A) implying that can also be directly ob¬ 
tained from 


dm 

dp 




(5) 


A similar scaling law is expected for r. The inset of 
Fig|^ shows the moment ratios around the critical point 
where the slope clearly increase (in absolute values) with 
size. The main plot shows the derivatives against size. 
Using the three methods, we estimate a critical exponent 
l/zz_L = 1.252(10), which is remarkably close to the expo¬ 
nent for Manna class l/zzx = 1.250(18) [2] and definitely 
ruling out the DP value l/zzx = 1.364(10) [T]. 


IV. CONCLUSIONS 

We investigate the effects of quenched coordination dis¬ 
order in continuous and discontinuous absorbing state 
phase transitions. Our extensive simulations of the ZGB 
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FIG. 6. (Color online) Determination of critical exponent 
for Manna model on VD lattices using different quantities 
X = \np, Inr and m. Inset: Moment ratio against control 
parameter around the critical point for L = 256, 512, 1024, 
and 2048. 

model on the VD lattice reveal the discontinuous nature 
of the absorbing state transition featured by the model 
remains unchanged under such a kind of disorder. Re¬ 
cently, it was shown that the Imry-Ma argument can be 
extended to non-equilibrium situations including absorb¬ 
ing states, and in addition, it was conjectured that first- 
order phase transitions cannot appear in low-dimensional 
disordered systems with an absorbing state. We showed 
that this is not always true: Our results for the ZGB 
model raise the interesting possibility that disorder in 
nonequilibrium APT may, under certain conditions, be 
irrelevant for the phase coexistence. The underlying rea¬ 
son for this is that the fluctuations induce correlated co¬ 
ordination disorder exhibited by the VD Lattice decay 
faster and are not able to preclude phase coexistence. 

In the case of continuous APT, we performed large- 
scale simulations of the A and Manna models on a 
Voronoi-Delaunay random lattice. Our results conhrm, 
as expected, that this kind of disorder does not alter 
the universality of the continuous transitions, support¬ 


ing that strong anticorrelations present in the VD ran¬ 
dom lattice makes topological disorder less relevant than 
uncorrelated randomness. 

Our findings corroborate a recent work of Barghatti 
and Vojta [IH] which shows systematically that the disor¬ 
der fluctuations of the VD lattice are featured by strong 
anticorrelations and decay faster than those of random 
uncorrelated disorder. In particular, it was shown that 
the random VD lattice has wandering exponent w = 
1/4 [48]. Hence, in this case, the Harris-Luck criterion 
yields that random connectivity is irrelevant at a clean 
critical point for v^_ > 2/3 that is satished for both 
Manna and DP universality classes. It is important to 
mention that in contrast to the A model, which belongs 
to DP class, even the strong disorder of uncorrelated lat¬ 
tice dilution (below the lattice percolation threshold) was 
found to be irrelevant for Manna class |32l|!^ . Therefore, 
our results are consistent with these findings, since the 
coordination disorder of VD lattice is weaker than lat¬ 
tice dilution. In addition, we determined the exponent 
l/v± = 1.252(10) for Manna class on VD lattice defi¬ 
nitely ruling out the DP value l/i^j. = 1.364(10). 

Further work should include the study of absorbing 
phase transitions on a three-dimensional random VD lat¬ 
tice, since it does not belong to the class of lattices with 
constrained total coordination [54]. In particular, accord¬ 
ing to the Harris criterion, the disorder might be relevant 
for the Manna class at least in three dimensions and there 
might be a dimensional difference between two and three 
dimensions. It would also be interesting to investigate if 
other kinds of correlated disorder are irrelevant for phase 
coexistence. 
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